#include "cuTypes.h"

namespace PerlinNoise
{
    SDH_CALL float  BrownianMotion   ( float3 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0, 
                                       float  gain       = 0.5 );

    SDH_CALL float  BrownianMotion   ( float4 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0, 
                                       float  gain       = 0.5 );

    SDH_CALL float3 BrownianMotionf3 ( float3 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0f, 
                                       float  gain       = 0.5f);

    SDH_CALL float3 BrownianMotionf3 ( float4 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0f, 
                                       float  gain       = 0.5f);
    SDH_CALL float  Turbulence       ( float3 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0f, 
                                       float  gain       = 0.5f);

    SDH_CALL float Turbulence        ( float4 point, 
                                       int    octaves    = 4, 
                                       float  lacunarity = 2.0f, 
                                       float  gain       = 0.5f);

    SDH_CALL float Random(int x, int y, int z, int t)
    {
	    int n = x + y * 59 + z * 263 + t * 12497;
	    n = (n<<13) ^ n;

	    return ( 1.0 - ( (n * (n * n * 3015 + 734368) + 123456789) & 0x7fffffff) / 1073741824.0f);
    }

    SDH_CALL float InterpolateCos(float x, float y, float a)
    {
	    float val = (1 - cos(a * 3.141592f)) * 0.5f;
	    return x * (1 - val) + y * val;
    }

    SDH_CALL float InterpolateS5(float x, float y, float a)
    {
        float val = a * a * a * (a * (a * 6.0f - 15.0f) + 10.0f);
        return x * (1 - val) + y * val;
    }

    SDH_CALL float InterpolateS3(float x, float y, float a)
    {
        float val = a * a * (3.0f - 2.0f * a) ;
        return x * (1 - val) + y * val;
    }

    SDH_CALL float Smooth2(float x, float y, int z, int t)
    {
	    float n1 = Random((int)x,     (int)y,     z, t);
	    float n2 = Random((int)x + 1, (int)y,     z, t);
	    float n3 = Random((int)x,     (int)y + 1, z, t);
	    float n4 = Random((int)x + 1, (int)y + 1, z, t);

	    float i1 = InterpolateS5(n1, n2, x - (int)x);
	    float i2 = InterpolateS5(n3, n4, x - (int)x);

	    return InterpolateS5(i1, i2, y - (int)y);
    }

    SDH_CALL float Smooth3(float x, float y, float z, int t)
    {
	    float i1 = Smooth2(x, y, (int)z,     t);
	    float i2 = Smooth2(x, y, (int)z + 1, t);

	    return InterpolateS5(i1, i2, z - (int)z);
    }

    SDH_CALL float Smooth4(float x, float y, float z, float t)
    {
	    float i1 = Smooth3(x, y, z, (int)t);
	    float i2 = Smooth3(x, y, z, (int)t + 1);

	    return InterpolateS5(i1, i2, t - (int)t);
    }

    SDH_CALL float Noise(float x, float y, float z)
    {
        return Smooth3(x, y, z, 0);
    }

    SDH_CALL float Noise (float3 p)
    {
        return Smooth3(p.x, p.y, p.z, 0);
    }

    SDH_CALL float Noise(float x, float y, float z, float w)
    {
        return Smooth4(x, y, z, w);
    }

    SDH_CALL float Noise(float4 p)
    {
        return Smooth4(p.x, p.y, p.z, p.w);
    }

    SDH_CALL float BrownianMotion(int x, int y, int z, int t)
    {
	    float noiseSum = 0.0f;

	    float freq = 0.015f;                                  // USER ADJUSTABLE
	    float gain = 0.65f;                                   // USER ADJUSTABLE
	    float octaves = 8;                                       // USER ADJUSTABLE
	    float amp   = 1;                                       // USER ADJUSTABLE

	    for(int lcv = 0; lcv < octaves; lcv++)
	    {
		    noiseSum = noiseSum + Noise(x * freq, y * freq, z * freq, t * freq) * amp;
		    freq = freq * 2;
		    amp = amp * gain;
	    }

	    float cloudCoverage = 0;                                     // USER ADJUSTABLE
	    float cloudDensity  = 1;                                     // USER ADJUSTABLE

	    noiseSum = (noiseSum + cloudCoverage) * cloudDensity;

	    if(noiseSum < 0) noiseSum = 0.0f;
	    if(noiseSum > 1) noiseSum = 1.0f;

	    return noiseSum;
    }


    SDH_CALL float BrownianMotion( float3 p, 
                                   int    octaves, 
                                   float  lacunarity, 
                                   float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;  

        for(int lcv=0; lcv<octaves; lcv++) 
        {
            noiseSum += Noise(p * freq) * amp;
            freq *= lacunarity;
            amp  *= gain;
        }

        return noiseSum;
    }

    SDH_CALL float BrownianMotion( float4 p, 
                                   int    octaves, 
                                   float  lacunarity, 
                                   float  gain)
    {  
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f; 

        for(int lcv=0; lcv<octaves; lcv++) 
        {
            noiseSum += Noise(p * freq) * amp;

            freq *= lacunarity;
            amp  *= gain;
        }

        return noiseSum;
    }

    SDH_CALL float3 BrownianMotionf3 ( float4 point1, 
                                       int    octaves, 
                                       float  lacunarity, 
                                       float  gain)
    {
        float3 f3Bm;

        float4 point2 = { 31, 78,  7, 3.7f};
        float4 point3 = {123, 11, 96, 5.1f};

        point2 = point2 + point1;
        point3 = point3 + point1;

        f3Bm.x = BrownianMotion(point1, octaves, lacunarity, gain);
        f3Bm.x = BrownianMotion(point2, octaves, lacunarity, gain);
        f3Bm.x = BrownianMotion(point3, octaves, lacunarity, gain);

        return f3Bm;
    }

    SDH_CALL float Turbulence( float3 p, 
                               int    octaves, 
                               float  lacunarity, 
                               float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;

        for(int lcv=0; lcv<octaves; lcv++) 
        {
            noiseSum += abs(Noise(p * freq)) * amp;

            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }

    SDH_CALL float Turbulence( float4 p, 
                               int    octaves, 
                               float  lacunarity, 
                               float  gain)
    {
        float freq = 1.0f, amp = 0.5f;
        float noiseSum = 0.0f;

        for(int lcv=0; lcv<octaves; lcv++) 
        {
            noiseSum += abs(Noise(p * freq)) * amp;

            freq *= lacunarity;
            amp  *= gain;
        }
        return noiseSum;
    }
}